#ifndef _HF_SOLVER_HO_H_
#define _HF_SOLVER_HO_H_

#include<iomanip>
#include<algorithm>//for sort algorithm
#include<cmath>
#include<eigen3/Eigen/Eigen>//matrix library

#include"orbital.h"
#include"system_ho.h"
#include"help.h"

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

using Eigen::ComplexEigenSolver;
using std::cout;
using std::stable_sort;
using std::abs;
using std::setw;



///spherical HF solver,for System_Tz_HO
class SHFSolver_HO
{
 public:
  typedef System_Tz_HO System_TzType;
  typedef typename System_TzType::DataType DataType;
  typedef typename System_TzType::OrbitalType OrbitalType;
  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,Dynamic,1> Vec;

 SHFSolver_HO():beta(0){}  
 SHFSolver_HO(const System_TzType * Ptr):pSystem_Tz(Ptr),beta(0){}
 SHFSolver_HO(const System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),A(_A),Z(_Z),realA(_A),beta(0){}


  void setbeta(double _beta)
  {
    beta=_beta;
  }
  
  void setup(const System_TzType * Ptr,int _A,int _Z)
  {
    pSystem_Tz=Ptr;
    A=_A;
    Z=_Z;
    realA=_A;
  }
  void setAZ(int _A,int _Z,int _realA=-1)
  {
    if(_realA==-1) realA=_A;
    else realA=_realA;
    A=_A;
    Z=_Z;
  }

  ///return onebody part under HO. basis of the hamiltonian
  DataType getOneBody(int bra,int ket) const
  {
    //return pSystem_Tz->get1B(realA,bra,ket);
    return pSystem_Tz->get1B(realA,bra,ket)*(1.-beta)+ beta*pSystem_Tz->get1B_p2_plus_r2(bra,ket);
  }
  ///return twobody part under HO. basis of the hamiltonian <ab->J|V|cd->J>,
  /*noting |ab->J> means sqrt(1+delta{ab})|ab,J>*/
  DataType getTwoBody(int a,int b,int c,int d,int J) const
  {
    //    return pSystem_Tz->get2B(realA,a,b,c,d,J);
    System_TzType::TwoBodyMatEle MatEle;
    
    pSystem_Tz->get2BmatAt(a,b,c,d,J,MatEle);
    if(MatEle.empty()) return 0;
    const double sqrt2=sqrt(2.);
    MatEle *= ( (a==b)? sqrt2 : 1.0 );
    MatEle *= ( (c==d)? sqrt2 : 1.0 );    
    return MatEle[0]+(beta-1)*MatEle[3]*pSystem_Tz->hbar_omega/realA-beta*MatEle[2]*pSystem_Tz->hbar_omega/realA;
  }
  
  ///1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole
  DataType SHFDiagram(int bra,int ket) const;

  ///return interaction part under HO. basis of the hamiltonian <ab->J|V|cd->J>,not include p1 \dot p2
  DataType getV(int a,int b,int c,int d,int J) const
  {
    return pSystem_Tz->getV(a,b,c,d,J);
  }
  ///return interaction mean field after HF. not include p1 \dot p2
  DataType getVHF(int bra,int ket) const;

  
  ///single particle hamiltonian, <bra|V1B|ket> + 1/(2j+1) sum_h{ (2J+1) < bra h->J | V2B | ket h->J >}, h is hole
  ///(group,a) corresponding to bra,(group,b) corresponding to ket
  DataType HFhamil(int group,int a,int b) const
  {
    int bra=pSystem_Tz->Groups[group][a];
    int ket=pSystem_Tz->Groups[group][b];
    return getOneBody(bra,ket)+SHFDiagram(bra,ket);
  }

  
  void printOrbitals(int num=-1) const
  {
    if(num==-1)
      num=HFOrbitals.size();
    cout<<"num"<<"\t"<<"l"<<"\t"<<"jj"<<"\t"<<"tz"<<"\t"<<"e"<<"\t"<<"OccuP"<<endl;    
    for(int i=0;i<num;i++)
      {
      	cout<<i<<"\t"<<HFOrbitals[i].l<<"\t"<<HFOrbitals[i].jj<<"\t"<<HFOrbitals[i].tz<<"\t"<<HFOrbitals[i].e<<"\t"<<HFOrbitals[i].OccuP<<endl;
      }
  }
  
  ///initialize
  void init();

  ///sort and fill orbitals.
  void fillOrbitals();

  ///construct matrix and diagonalize, generate new HF orbitals
  void diag();
  
  ///iteration
  void iter();

  ///cal. HF energy
  void calEnergy();



  double beta;
  
  int A,Z;//core 
  int realA;//the A in Hamiltonian
  vector< SHF_Orbital<DataType> > HFOrbitals;
  DataType Energy;

  const System_TzType * pSystem_Tz;
};

//*********************************************************************
#endif
